library(factoextra)
library(pheatmap)
library(tidyverse)
theme_set(theme_bw())

We usually think of high-dimensional data as consisting of multiple measures on a group of samples:

• Number of “reads” of different bacterial proteins from a set of soil samples
• Decathlon scores from different competitors
• Health measures from different children

Many scientists traditionally think of high-dimensional data as having parallel, continuous measures:

• read matrices from soil samples

These may be complemented by a smaller number of “metadata” variables, which may be more diverse in type (count, categorical, etc.):

• environmental variables associated with soil samples

More and more datasets don’t follow this:

• Canadian longitudinal study on aging has a huge number of variables per person with a wide mixture of types

We study the rows (samples) using the columns (measures)

• What do the observed proteins tell us about the functional relationships between different soil samples?
• What does differential success in decathlon events tell us about the athletes?

But we can also do the opposite!

• What do measurements across soil samples tell us about the functional relationships between proteins?
• What does differential success of athletes tell us about the relationship between events?

print(names(decathlon2))  ## from factoextra pkg

##  [1] "X100m"        "Long.jump"    "Shot.put"     "High.jump"   
##  [5] "X400m"        "X110m.hurdle" "Discus"       "Pole.vault"  
##  [9] "Javeline"     "X1500m"       "Rank"         "Points"      
## [13] "Competition"

dec_frame <- (decathlon2[1:10]
        %>% as_tibble()
    %>% rename_all(sub, pattern="^X", replacement="Run_")
    %>% mutate_at(vars(contains('Run_')), ~(-1*.))
    %>% rename_all(sub, pattern="110m.", replacement="")
    %>% mutate_all(~c(drop(scale(.))))
)

Iris pairs plot (iris.R)

• Slow, and not very good for more than 5 or 6 variables
• Base R may be faster for quick viz

dec_mat <- as.matrix(dec_frame)
heatmap(dec_mat)

heatmap(t(dec_mat))

pheatmap(cor(dec_mat), cell.width = 10, cell.height = 10)

Better for visualizing groups of events

pheatmap(cor(t(dec_mat)), cell.width = 4, cell.height = 4)

Better for visualizing groups of events

A beautiful decomposition based on the idea that data points are points in a Euclidean space

• Need to think about scaling

We can think about the PCA as a decomposition (making observed points from idealized points)

• And relax it by requiring non-negative combinations of non-negative components (NMF)

Or we can think about it as minimizing distances:

• And relax it with non-Euclidean distances (PCoA)
• … or a rank-based approach (NMDS)

pca_ath <- prcomp(dec_mat, scale=TRUE)
fviz_screeplot(pca_ath)

fviz_pca_ind(pca_ath)

fviz_pca_var(pca_ath)

fviz_pca_var(pca_ath, axes=c(2, 3))

View scores and components:

fviz_pca_biplot(pca_ath)

load <- pca_ath$sdev * pca_ath$rotation
pheatmap(load, cluster_cols=FALSE)

WRONG!!!!

load <- pca_ath$sdev * pca_ath$rotation
pheatmap(load, cluster_cols=FALSE)

WRONG!!!!

library(vegan)

## Loading required package: permute

## Loading required package: lattice

## This is vegan 2.5-6

mds <- metaMDS(dec_mat, distance="euclidean")

## 'comm' has negative data: 'autotransform', 'noshare' and 'wascores' set to FALSE

## Run 0 stress 0.175356 
<<<<<<< HEAD
## Run 1 stress 0.175356 
## ... Procrustes: rmse 3.75291e-06  max resid 1.004819e-05 
## ... Similar to previous best
## Run 2 stress 0.2029784 
## Run 3 stress 0.175356 
## ... Procrustes: rmse 1.355825e-05  max resid 3.876792e-05 
## ... Similar to previous best
## Run 4 stress 0.2028176 
## Run 5 stress 0.175356 
## ... Procrustes: rmse 4.319605e-05  max resid 0.0001589884 
## ... Similar to previous best
## Run 6 stress 0.1975353 
## Run 7 stress 0.2035142 
## Run 8 stress 0.1787377 
## Run 9 stress 0.2020938 
## Run 10 stress 0.1787376 
## Run 11 stress 0.2146083 
## Run 12 stress 0.2035077 
## Run 13 stress 0.1783895 
## Run 14 stress 0.1975374 
## Run 15 stress 0.1754619 
## ... Procrustes: rmse 0.02261613  max resid 0.08184932 
## Run 16 stress 0.1950854 
## Run 17 stress 0.2028353 
## Run 18 stress 0.1950854 
## Run 19 stress 0.1788135 
## Run 20 stress 0.1787378 
## *** Solution reached

plot(mds)

## species scores not available

The biplot is hard (but we can help you if you want to do it)

library(GGally)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

## 
## Attaching package: 'GGally'

## The following object is masked from 'package:dplyr':
## 
##     nasa

The biplot is hard (but we can help you if you want to do it)

library(GGally)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

## 
## Attaching package: 'GGally'

## The following object is masked from 'package:dplyr':
## 
##     nasa

ggparcoord(dec_frame)

library(directlabels)
ath_frame <- (
    as.data.frame(t(as.matrix(dec_frame)))
    %>% mutate(event=names(dec_frame))
)
ggparcoord(ath_frame)

ggparcoord(ath_frame, groupColumn="event")

library(andrews)
<<<<<<< HEAD
andrews(dec_frame)

andrews(ath_frame)

andrews(ath_frame)

Not implemented; may suffer similar problems

Many of these charts might work if the variables were in a suitable order:

<<<<<<< HEAD
• heatmap provides an order =======

• heatmap provides an order

  >>>>>>> refs/remotes/origin/gh-pages
  • hierarchical

• other ordering methods take a more holistic view

<<<<<<< HEAD
• PCA and other decompositions could also provide an order! =======

• PCA and other decompositions could also provide an order!

  >>>>>>> refs/remotes/origin/gh-pages
  • Should we combine PCA with Andrews or parcoord?